# !/usr/bin/env python
# -*- coding: utf-8 -*-
"""
@Time        : 2021/5/8 15:09
@Author      : Albert Darren
@Contact     : 2563491540@qq.com
@File        : RGB_split.py
@Version     : Version 1.0.0
@Description : TODO
@Created By  : PyCharm
"""
import numpy as np

# 自己原创平方根法
def cholesky_decomposition(sym_pos_define_matrix: np.ndarray, right_hand_side_vector: np.ndarray):
    rows, columns = sym_pos_define_matrix.shape
    # judge if it is a square matrix
    if rows == columns:
        # 判断是否对称
        if (sym_pos_define_matrix.T == sym_pos_define_matrix).all():
            for k in range(rows):  # 判断各阶顺序主子式是否为0,即判定是否正定
                if np.linalg.det(sym_pos_define_matrix[:k + 1, :k + 1]) <= 0:
                    raise Exception("cannot decompose")
            else:
                # 初始化单位下三角矩阵L
                square_ones_matrix = np.ones((rows, columns))
                lower_triangle_matrix = np.tril(square_ones_matrix)
                for j in range(rows):
                    prod = 0
                    for k in range(j):
                        prod += lower_triangle_matrix[j, k] * lower_triangle_matrix[j, k]
                    lower_triangle_matrix[j, j] = np.sqrt(sym_pos_define_matrix[j, j] - prod)

                    for i in range(j + 1, rows):
                        prod = 0
                        for k in range(j):
                            prod += lower_triangle_matrix[i, k] * lower_triangle_matrix[j, k]
                        lower_triangle_matrix[i, j] = (sym_pos_define_matrix[i, j] - prod) / lower_triangle_matrix[j, j]

        else:
            raise Exception("error,the input matrix must be a symmetric-matrix")
    else:
        raise Exception("ERROR:please pass a square matrix.")
    return np.linalg.inv(lower_triangle_matrix.transpose()) @ np.linalg.inv(
        lower_triangle_matrix) @ right_hand_side_vector


if __name__ == '__main__':
    import cProfile
    symmetric_positive_define_matrix2 = np.array([[4, 2, -4, 0, 2, 4, 0, 0],
                                                  [2, 2, -1, -2, 1, 3, 2, 0],
                                                  [-4, -1, 14, 1, -8, -3, 5, 6],
                                                  [0, -2, 1, 6, -1, -4, -3, 3],
                                                  [2, 1, -8, -1, 22, 4, -10, -3],
                                                  [4, 3, -3, -4, 4, 11, 1, -4],
                                                  [0, 2, 5, -3, -10, 1, 14, 2],
                                                  [0, 0, 6, 3, -3, -4, 2, 19]], dtype=np.float64)
    column_vector = np.array([0, -6, 20, 23, 9, -22, -15, 45], dtype=np.float64).reshape((8, 1))
    # print(cholesky_decomposition(symmetric_positive_define_matrix2, column_vector))
    cProfile.run("cholesky_decomposition(symmetric_positive_define_matrix2, column_vector)",
                 filename="result.out", sort="cumulative")

